function Plot_Entropy(parms,gesol,stats,prob_Nxz)

    m = 2; n = 3; idxplot = 0;

    idx1y = find(gesol.entropy.horizon==12);
    idx5y = find(gesol.entropy.horizon==60);
    
    % should only be plotted for regions that's actually visited
    prob_N = sum(sum(prob_Nxz,3),2);
    cdf_N = cumsum(prob_N);
    idx_N1 = find(cdf_N>=0.0001,1,'first');
    idx_N2 = find(cdf_N<=0.9999,1,'last');
    
    prob_x = sum(sum(prob_Nxz,3),1);
    cdf_x = cumsum(prob_x);
    idx_x1 = find(cdf_x>=0.0001,1,'first');
    idx_x2 = find(cdf_x<=0.9999,1,'last');
    
    IDX_N_SELECT = idx_N1:idx_N2; IDX_X_SELECT = idx_x1:idx_x2;
    phi_of_x = @(x)1./(1+exp(-x));
    [NN1,PP1] = meshgrid(parms.grid_N(IDX_N_SELECT),phi_of_x(parms.grid_x(IDX_X_SELECT)));
    
    figure
    
    idxplot = idxplot + 1; subplot(m,n,idxplot)
    mesh(NN1,PP1,gesol.entropy.CondEntropy(IDX_N_SELECT,IDX_X_SELECT,1,idx1y)')
    xlabel('N')
    ylabel('\phi')
    title('CondEntropy per year, 1 year, z=1')
    axis square
    
    idxplot = idxplot + 1; subplot(m,n,idxplot)
    mesh(NN1,PP1,gesol.entropy.CondEntropy(IDX_N_SELECT,IDX_X_SELECT,1,idx5y)'/5)
    xlabel('N')
    ylabel('\phi')
    title('CondEntropy per year, 5 year, z=1')
    axis square
    
    idxplot = idxplot + 1; subplot(m,n,idxplot)
    mesh(NN1,PP1,gesol.entropy.CondEntropy(IDX_N_SELECT,IDX_X_SELECT,1,idx1y)'...
        -gesol.entropy.CondEntropy(IDX_N_SELECT,IDX_X_SELECT,1,idx5y)'/5)
    xlabel('N')
    ylabel('\phi')
    title('1 year - 5 year, z=1')
    axis square
    
    idxplot = idxplot + 1; subplot(m,n,idxplot)
    mesh(NN1,PP1,gesol.entropy.CondEntropy(IDX_N_SELECT,IDX_X_SELECT,2,idx1y)')
    xlabel('N')
    ylabel('\phi')
    title('CondEntropy per year, 1 year, z=2')
    axis square
    
    idxplot = idxplot + 1; subplot(m,n,idxplot)
    mesh(NN1,PP1,gesol.entropy.CondEntropy(IDX_N_SELECT,IDX_X_SELECT,2,idx5y)'/5)
    xlabel('N')
    ylabel('\phi')
    title('CondEntropy per year, 5 year, z=2')
    axis square
    
    idxplot = idxplot + 1; subplot(m,n,idxplot)
    mesh(NN1,PP1,gesol.entropy.CondEntropy(IDX_N_SELECT,IDX_X_SELECT,2,idx1y)'...
        -gesol.entropy.CondEntropy(IDX_N_SELECT,IDX_X_SELECT,2,idx5y)'/5)
    xlabel('N')
    ylabel('\phi')
    title('1 year - 5 year, z=2')
    axis square
    
    figure
    
    idxplot = 0;
    
    idxplot = idxplot + 1; subplot(m,n,idxplot)
    mesh(NN1,PP1,gesol.entropy.E_dlogM(IDX_N_SELECT,IDX_X_SELECT,1,idx1y)')
    xlabel('N')
    ylabel('\phi')
    title('E_t[m(t+T)-m(t)] per year, 1 year, z=1')
    axis square
    
    idxplot = idxplot + 1; subplot(m,n,idxplot)
    mesh(NN1,PP1,gesol.entropy.E_dlogM(IDX_N_SELECT,IDX_X_SELECT,1,idx5y)'/5)
    xlabel('N')
    ylabel('\phi')
    title('E_t[m(t+T)-m(t)] per year, 5 year, z=1')
    axis square
    
    idxplot = idxplot + 1; subplot(m,n,idxplot)
    mesh(NN1,PP1,gesol.entropy.E_dlogM(IDX_N_SELECT,IDX_X_SELECT,1,idx1y)'...
        -gesol.entropy.E_dlogM(IDX_N_SELECT,IDX_X_SELECT,1,idx5y)'/5)
    xlabel('N')
    ylabel('\phi')
    title('1 year - 5 year, z=1')
    axis square
    
    idxplot = idxplot + 1; subplot(m,n,idxplot)
    mesh(NN1,PP1,gesol.entropy.E_dlogM(IDX_N_SELECT,IDX_X_SELECT,2,idx1y)')
    xlabel('N')
    ylabel('\phi')
    title('E_t[m(t+T)-m(t)] per year, 1 year, z=2')
    axis square
    
    idxplot = idxplot + 1; subplot(m,n,idxplot)
    mesh(NN1,PP1,gesol.entropy.E_dlogM(IDX_N_SELECT,IDX_X_SELECT,2,idx5y)'/5)
    xlabel('N')
    ylabel('\phi')
    title('E_t[m(t+T)-m(t)] per year, 5 year, z=2')
    axis square
    
    idxplot = idxplot + 1; subplot(m,n,idxplot)
    mesh(NN1,PP1,gesol.entropy.E_dlogM(IDX_N_SELECT,IDX_X_SELECT,2,idx1y)'...
        -gesol.entropy.CondEntropy(IDX_N_SELECT,IDX_X_SELECT,2,idx5y)'/5)
    xlabel('N')
    ylabel('\phi')
    title('1 year - 5 year, z=2')
    axis square
    
    figure
    
    m = 2; n = 6; idxplot = 0;
    
    idxplot = idxplot + 1; subplot(m,n,1:2)
    hold on
    plot(stats.entropy.horizon(:)/12,stats.entropy.MeanCondEntropy(:)...
        ./(stats.entropy.horizon(:)/12),'-');
    xlabel('horizon, years')
    title('E[CondEntropy]/Horizon')
    ylabel('Average entropy per year')
    axis square
    
    idxplot = idxplot + 1; subplot(m,n,3:4)
    hold on
    plot(stats.entropy.horizon(:)/12,stats.entropy.MeanCondEntropy(:)...
        ./(stats.entropy.horizon(:)/12),'-');
    plot(stats.entropy.horizon(:)/12,stats.entropy.MeanCondEntropy_given_z...
        ./(stats.entropy.horizon(:)/12),'-.');
    xlabel('horizon, years')
    title('E[CondEntropy]/Horizon')
    ylabel('Average entropy per year')
    axis square
    lgdtxt = cell(parms.Nz+1,1);
    lgdtxt{1} = 'unconditional';
    for iz = 1:parms.Nz
        lgdtxt{1+iz} = ['conditional: z=',num2str(iz)];
    end
    legend(lgdtxt)
    legend boxoff
    
    idxplot = idxplot + 1; subplot(m,n,5:6)
    hold on
    plot(stats.entropy.horizon(:)/12,stats.entropy.EdlogM_given_z...
        ./(stats.entropy.horizon(:)/12),'-');
    xlabel('horizon, years')
    title('E[m(t+H)-m(t)|z(t)]/Horizon')
    ylabel('Average entropy per year')
    axis square
    lgdtxt = cell(parms.Nz,1);
    for iz = 1:parms.Nz
        lgdtxt{iz} = ['conditional: z=',num2str(iz)];
    end
    legend(lgdtxt)
    legend boxoff
    
    idxplot = idxplot + 1; subplot(m,n,7:9)
    [TT2,NN2] = meshgrid(stats.entropy.horizon(:)/12, parms.grid_N(IDX_N_SELECT));
    mesh(TT2,NN2,(stats.entropy.MeanCondEntropy_given_N(:,IDX_N_SELECT)./(repmat(stats.entropy.horizon(:),[1 numel(IDX_N_SELECT)])/12))')
    xlabel('horizon, years')
    title('E[CondEntropy|N]/Horizon')
    ylabel('N')
    
    idxplot = idxplot + 1; subplot(m,n,10:12)
    [TT3,PP3] = meshgrid(stats.entropy.horizon(:)/12, phi_of_x(parms.grid_x(IDX_X_SELECT)));
    mesh(TT3,PP3,(stats.entropy.MeanCondEntropy_given_x(:,IDX_X_SELECT)./(repmat(stats.entropy.horizon(:),[1 numel(IDX_X_SELECT)])/12))')
    xlabel('horizon, years')
    title('E[CondEntropy|\phi]/Horizon')
    ylabel('\phi')
    
    figure
    
    m = 1; n = 2; idxplot = 0;
    
    idxplot = idxplot + 1; subplot(m,n,idxplot)
    hold on
    NT = stats.entropy.horizon(end);
    temp = permute(repmat(reshape(stats.entropy.MeanCondEntropy_given_z(12,:),[],1),[1 NT]),[2 1]) ...
        - stats.entropy.MeanCondEntropy_given_z...
        ./(stats.entropy.horizon(:)/12);
    plot(stats.entropy.horizon(:)/12,temp,'-');
    xlabel('horizon, years')
    title('Slope component 1: L_t^{(1)} - L_t^{(n)}/n')
    axis square
    lgdtxt = cell(parms.Nz,1);
    for iz = 1:parms.Nz
        lgdtxt{iz} = ['conditional: z=',num2str(iz)];
    end
    legend(lgdtxt)
    legend boxoff
    grid on
    set(gca,'xlim',[1 5])
    
    idxplot = idxplot + 1; subplot(m,n,idxplot)
    hold on
    NT = stats.entropy.horizon(end);
    temp = permute(repmat(reshape(stats.entropy.EdlogM_given_z(12,:),[],1),[1 NT]),[2 1]) ...
        - stats.entropy.EdlogM_given_z...
        ./(stats.entropy.horizon(:)/12);
    plot(stats.entropy.horizon(:)/12,temp,'-');
    xlabel('horizon, years')
    title('Slope component 2: E_t[m(t+1)-m(t)] - E_t[m(t+n)-m(t)]/n')
    axis square
    lgdtxt = cell(parms.Nz,1);
    for iz = 1:parms.Nz
        lgdtxt{iz} = ['conditional: z=',num2str(iz)];
    end
    legend(lgdtxt)
    legend boxoff
    grid on
    set(gca,'xlim',[1 5])

end